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y '. Abstract 

-^ , Using network models consisting of gap junction coupled Wang-Buszaki neurons, we demonstrate 

r^ ' that it is possible to obtain not only synchronous activity between neurons but also a variety of 



constant phase shifts between and tt. We call these phase shifts intermediate stable phase- 



m , locked states. These phase shifts can produce a large variety of wave-like activity patterns in 

> ■ 
^C I one-dimensional chains and two-dimensional arrays of neurons, which can be studied by reducing 



the system of equations to a phase model. The 27r periodic coupling functions of these models are 



m . characterized by prominent higher order terms in their Fourier expansion, which can be varied by 

O 

^N I changing model parameters. We study how the relative contribution of the odd and even terms 

affect what solutions are possible, the basin of attraction of those solutions and their stability. 

K> ■ These models may be applicable to the spinal central pattern generators of the dogfish and also to 

- - - the developing neocortex of the neonatal rat. 



I. INTRODUCTION 



In the past few decades, neuroscientists have discovered the importance of oscillations in 
the brain. Oscillations occur in many neural networks and play key roles in sensory as well 



as motor systems 



ll 



y. On 



the smallest scale of the nervous system, individual neurons can 



behave as complex nonlinear oscillators. [3| On a larger scale, oscillatory patterns in EEG 
(electroencephalographic) recordings have been shown to correspond to different cognitive 
statesjSl]. Phase oscillator models (obtained from systems of weakly coupled oscillators) have 
proven extremely useful in the analysis of such systems [^ |5| . 

A pair of symmetric weakly coupled oscillators always has at least two possible periodic 
(phase-locked) states: the synchronous state and the anti-phase atate (where they are a half 
a cycle out of phase). Suppose that for some range of parameters the synchronous solution 
is stable and as one of these parameters is varied, the synchronous solution loses stability. 
In a symmetrically coupled system the synchronous solution generically loses stability in a 
pitchfork bifurcation. When this pitchfork bifurcation is supercritical it gives rise to two 
new stable phase-locked solutions in which the phase difference between oscillators is some 
constant between and vr. We refer to such solutions as Intermediate stable phase locked 
states. Intermediate stable phase locked states between pairs of neurons can be produced 
in a variety of neuron models by using a variety of coupling schemes as well as adjusting 



the parameters which affect the excitability of individual neurons |6|][7|[8[. The first work 
examining this behavior in conductance-based models was by Vanvreeswijk et al [8]. In this 
paper, the authors studied integrate-and-fire as well as Hodgkin-Huxley neurons. Cymbalyuk 
has both modeled and experimentally demonstrated an intermediate stable phase-locked 
state in a system of two "silicon neurons" 9| 10|]. These intermediate phase-locked states 
can lead to wave-like behavior in large networks. 

Wave activity is ubiquitous in the brain. Traveling waves of electrical activity in the 
brains of animals have been observed in a variety of species and occur in a diverse set 
of structures. These structures include the retina, olfactory cortex Oil, peri-geniculate 

n n n nnrnnn 

nucleus [121], neocortex |13| and spinal cord [l[, among others [1^ , [12[ , [I3j , [111 , [15[ [12| . In 
phase models, waves correspond to constant nonzero phase-differences between successive 
pairs of oscillators. There are many ways to generate such phase differences (see [l6[ , section 
8.3.5) such as a gradient in natural frequencies J], pacemakers, or manipulation of the 



boundary conditions 17|]. In some swimming organisms, it is very important that a fixed 



phase-lag be maintained over a wide range of frequencies. For example, in the lamprey, the 



lag is about 27r/100 [? ], while for the crayfish it is 271 /A [19j. As we will show in this paper, 
intermediate stable phase-locked states provide a simple way to produce waves with a stable 
fixed inter-oscillator phase- difference. 

The paper is organized as follows. We first show that by varying the excitability of gap- 
junction coupled Wang-Buszaki neurons, we can produce intermediate stable phase-locked 
states. We then explore the consequences of these states in one-dimensional nearest-neighbor 
coupled chains. We prove the existence and stability of a wide variety of complex waves 
including simple waves and zig-zag waves. We then study two-dimensional arrays and find 
two-dimensional analogues of zig-zag and regular waves. 



II. GAP JUNCTION COUPLING BETWEEN PAIRS OF WANG-BUSZAKI NEU- 
RONS 

A. Measuring the Phase Difference Between Spiking Neurons 

The Wang-Buszaki model is a conductance based neuron model derived from the Hodgkin 



Hux! 



ey 



model. It was originally used to describe fast spiking interneurons in the hippocam- 



pus 20|. In this section we examine pairs of Wang-Buszaki neurons reciprocally coupled with 
gap junctions and ask what happens as we vary parameters which affect the excitability of 
the individual cells. Gap junctions are specialized ion channels connecting the cytoplasm 
of the pre and post-synaptic cell. A depolarizing ionic current is driven by the potential 
difference between the cells. In order to understand the role that gap junctions play in 
rhythmically oscillating networks, consider just two neurons coupled with gap junctions. 
The coupled neuron equations will have the following form: 

CVl = -iNa -h-lL-h + 9{V2 - Vi) 
CV2 = -iNa -Ik-lL-l0 + 9{Vl - V2) 

(1) 



These currents are given by the equations: 

iNa = gNaml^h{V - Vncl) 

h = 9kn\V-Vk) 

h = gL{y-VL) (2) 

C is the the membrane capacitance measured in units of fiF/cm'^. The parameters: V^ayVk 
and Vl are the reversal potentials for the ion channels. The parameters g, qng, Qk and ql are 
constants describing the conductances of their respectively labeled ion channels. Typically, 
they are measured in units of mS/cm? . The m and the h are the time-dependent activation 
and decativation variables which are described by the equations: 

'^ = ^{a^(y){l - h) - p^{V)h) 
^ = vicymiV)il - m) - PUV)m) (3) 

In these equations, rj is an overall channel switching rate determined by the temperature 



of the system 211]. The nonlinearities ai{V),Pi{V), as well as the parameters used in the 
simulations are referenced in the Appendix El 

Previously, Pfeuty et al. have studied pairs of Wang-Buszaki neurons with gap junction 
coupling and shown that by varying the potassium and sodium conductances, one can vary 
the stability of the synchronous in-phase and anti-phase states [22]. In order to find pa- 
rameter regimes in which intermediate phase- locked states exist, we varied the potassium 
conductance, Qk and also the temperature dependent rate constant, r]. Parameters such as 
these play key roles in the behavior of both central pattern generators and large scale oscilla- 



tory networks 23|] 2J]. The most important consequence of varying these parameters is that 



the absolute refractory period decreases. The neuron is less hyperpolarized after the action 



potential for larger rj or smaller g^ |22(]. In other words, this smaller absolute refractory 
region allows for the neuron to fire more quickly [22]. Consider the upper half of Figure [TJ 
Panels A and B are calculations of the phase difference between two Wang-Buszaki neurons 
computed after 1500 ms of integration time for different values of Qk (the maximal potassium 
conductance) and rj (the temperature dependent time constant). The phase difference, (p 
was measured as the difference between the times at which Vi and V2 cross zero with a posi- 
tive slope divided by the period of the oscillation. Holding the stimulus current constant at 
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FIG. 1. A. The phase difference between two Wang Buszaki neurons (Fuh model) as a function of 
the parameter rj. The diagram clearly illustrates a pitchfork bifurcation connecting the synchronous 
and anti-phase solutions. B. Plot of the phase difference for varying values of gk- C. Plot of the 
odd part of the interaction function (cf equations (j8|9p ) for different values of rj (the dimensionless 
temperature dependent time constant). D. Calculation of the odd part of the interaction function 
for different values of g^. The zeros of these functions correspond to the stable phase locked 
solutions. 



io = 0.63nA and varying either the parameter t] ot gt, we are able to demonstrate a super- 
critical pitchfork bifurcation in the phase difference between neurons. Figure [TJ\ illustrates 
the bifurcation of the system of two neurons from synchrony to anti-phase behavior as t] is 
increased. For values of rj between 5.0 and 7.0 the system passes through all possible stable 
relative phase-locked solutions varying between 0.0 and vr. Figure [T]B is a plot of the phase 
difference between neurons as we decrease Qk- As we decrease gk between 8.0 and 5.0 we 
also obtain intermediate stable phase-locked solutions. 



B. Calculation Of The Interaction Function 

In order to simplify the analysis, we reduce our system of neurons to a phase model. 
This can be accomplished by applying Malkins Theorem [2|. Consider two weakly coupled 
systems (By weak, we mean_that the coupling acts to only affect the phase and not the 
amplitude of the oscillation [2j]) of coupled differential equations describing tonically firing 
neurons: 

V; = F{V,) + eCi(\/i, V2) + 0{e') (4) 

V^ = F{V2) + eC2{V2, Vi) + 0{e^). (5) 

Then Malkin's theorem states: 

Vi{t) = VoiOi) + 0(e), V2{t) = ¥0(62) + 0(e) 

where Vq is the solution to the homogenous equation. Furthermore: 

-^ = l + eH,ie2-9,) + Oie') (6) 

-^ = l + eH2{ei-e2) + 0{^) (7) 

where: 

HM--=^I V*it)-CAVoit),Voit + <l))]dt. (8) 

Here V* is the adjoint of the linearized equation (jl]). Let (p := 62 — Oi- Then: 

^ = -2eH,,M := e[/J(-0) - iJ(0)]. (9) 

The right hand sides of these phase model equations are described by interaction functions 
which are derived from the full model. Zeros of Hodd{<P) determine the possible phase-locked 
states and those for which H'^^{(j)) > are stable. 

We calculated the interaction functions for several values of gt and 77. Figure IC and 
ID shows the odd portions of the interaction functions calculated with various parameter 
values. We write these interaction functions in terms of their Fourier series: 

H{x) = — + 2_. \ '^n cos{nx) + bn sin(nx) J (10) 

For example. Table [T] shows the first few Fourier modes of the interaction function computed 
from the full model with rj = 6. This table demonstrates the substantial contribution of 
higher order Fourier terms to the interaction function near the bifurcation point. 



ao = 5.1974931 
ai = -2.9970722 h = 0.47408548 

02 = -0.92187762 62 = -0.36833799 

03 = -0.44113794 63 = -0.2577318 

04 = -0.25482759 64 = -0.15762125 
as = -0.16416954 65 = -0.09083201 
ae = -0.11295291 ftg = -0.048487604 

TABLE I. Fourier Coefficients of the Interaction function, equation ([8]) (computed from the full 
model with ry = 6). Near the bifurcation point from synchrony there are substantial higher order 
even and odd Fourier terms. 

Observe that the first two odd Fourier terms in Tabled are the dominant ones (61, 62)- As 
the shape of Hodd{<P) depends only on the odd terms, aj have no effect on the existence and 
stability of the locked solutions for a pair of symmetrically coupled oscillators. However, once 
there are more than just two oscillators, the even terms play an important role, particularly 
as fasr as stability is concerned. We find that as rj changes between 5 and 7, the Fourier 
coefficient, 62 remains negative and 61 goes to zero. For this reason, in the rest of this paper, 
we will use a truncated version of the H function that contains only three terms: 

H{x) = bi sin(x) + 62 sin(2x) + ai cos(a;) (11) 

Unless otherwise stated, we typically take 61 = 1 and 62 = —.75, giving zeros of Hodd at 
(f) = 0, vr, ±cos~^(2/3) ~ ±0.847. We remark that this particular form of Hoddi'P) arises in 
the bead on a hoop instability [25|(Chapt 3.5) and in a model for the coordination of finger 



tapping ((26j). If 62 < is fixed and negative, then as hi decreases from a large positive value, 
the synchronous solution loses stability (at 61 = —262) and a branch of intermediate stable 
phase-locked solutions bifurcates. This branch remains stable until hi becomes sufficiently 
negative hi < 262 whereupon the anti-phase solution is stable. 

III. WAVES IN LARGE NETWORKS 

This section is a study of wave behavior in both chains and two dimensional arrays of 
neurons with nearest neighbor coupling in regimes where there is an intermediate stable 



phase-locked state for pair-wise coupling. Primarily, we study phase models which use the 
interaction functions (or approximations of them) derived from the Wang-Buszaki model. 
Our models may be relevant to patterns of wave activity in the neonatal rat. Peinado et 
al. was able to observe wave activity in gap-junction coupled interneurons in rat neocortex 
prior to day twelve of development 13|]. Furthermore, he was able to enhance these waves by 
applying halothane and picrotoxin. Picrotoxin blocks inhibitory synapses while halothane 
reduces the potassium conductance. In general, he observed that the reduction in potassium 
directly led to the formation of waves. Since our intermediate stable phase-locked states can 
occur in gap junction coupled neurons by reducing the potassium conductance, this may be 
experimental evidence that this effect plays a role in wave formation in a two-dimensional 
network. 

We focus on a specific type of solution to the phase equations known as anti-waves. Anti- 
waves were first studied in two papers by Ermentrout and Kopell 17| 27||. Similar phenomena 
have been examined by Strogatz et al., (uniformly twisted waves) and Blasius et al., (ID 
quasiregular concentric waves) 28|] 29| 30|. Anti- waves consist of waves either initiated at 
the ends and colliding in the middle or waves initiating from the middle and terminating 
at the ends. The latter are, in a sense, equivalent to one- dimensional target patterns. The 
wavenumber for these waves (the spatial gradient of the phase) shows an abrupt change of 
sign which we will call a kink. Similar phase waves have been demonstrated in mechanical 



systems 



311]. Anti- waves have been experimentally observed in the spinal cords of dogfish 



32| and may well be present in other biological tissue. For instance, similar patterns of 



electrical activity have been observed in the muscle of the colon of a cat 33| . Central 



pattern 



generators in the fins of electric fish have also been known to produce anti- waves 271]. These 
animals are able to produce a variety of complex waves in which the "kink" or lead oscillator 
in the wave is able to shift |2 71]. 

We want to stress that our mechanism for generating anti- waves is fundamentally different 
than previous papers. Furthermore, our anti-waves can form anywhere along the chain 
and are much more robust to perturbations than in previous models 17|. The rest of this 
paper is organized into several sections. We begin in section IIIIAI by discussing the basic 
phase models and boundary conditions. This is followed by an analysis of both an ordinary 
traveling wave ( ]IIIB|) and anti- wave solutions (1111 C|) . In section IIIIDt we demonstrate that 
the probability of obtaining a particular solution depends on the relative contribution of 



8 



the even component of the interaction function. We show that starting from the anti-wave 
solution, if the even component is sufficiently large, perturbations initiated at one end of 
the chain can propagate down the chain and shift the position of the kink. Finally, inlVlwe 
demonstrate that this analysis can be extended to higher dimensions and that a variety of 
anti-wave patterns are possible in a two dimensional oscillator arrays. 



A. Models And Boundary Conditions 



The models we consider were introduced by Kopell and Ermentrout in a 1986 paper 17 1. 
These models primarily describe networks of neurons with nearest neighbor coupling. In 
analyzing these equations we apply two types of boundary conditions: periodic boundary 
conditions and non-refiecting boundary conditions. For a system of A^ + 1 neurons with 
periodic boundary conditions, the system of phase equations may be written: 



e'l = wi + HLieM+^ - 9,) + Hni92 - 9^) 
92 = UJ2 + Hl{9i - 92) + Hr{9:, - 92) 

9n+i = ^N+i + Hl{9n — 9]^^i) + Hji{9i — 6'jv^i). 

(12) 



In these equations, the Ui represents the natural frequencies of the oscillators. We denote 
the coupling in the two possible directions as Hl{(J)) and Hr{(J)). In general, we assume that 
the coupling is isotropic, so that: Hji{(j)) = Hl{(J)) = H{(j)). Ultimately, we are interested in 
phase-locking behavior, thus we make the change of variables: 0j = 0j+i — 0j. This results 
in a system of A^ phase equations: 

9 



N 
62 = /\UJ2 + ^(03) + ^(-02) - i^(-0l) - i^(02) 



N 



(13) 



In these equations Awj is the frequency gradient between oscillators. In most of our simu- 
lations, we assume identical frequencies so Awj = 

The second type of boundary condition that we use is a variation of what are known as non- 
reflecting boundary conditions |34| . Non- reflecting boundary conditions are implemented in 
order to attempt to eliminate reflections and to "trick" the neurons at the ends of the chains, 
neuron 1 and neuron A^ -|- 1, into behaving as though the chain is infinite. If one were to 
think of the chain as being a continuous system, then the boundary conditions are simply a 
statement that ^^' ' = when evaluated at the ends of the chain. There is a precedent for 



using such boundary conditions in nonlinear oscillator problems, for an example, see 311] • 
Applying these boundary conditions to our phase equations, we have[35|: 

9o = O2 
On+1 = On-1, (14) 

and our phase equations become: 

01 = Aui + H{<j)2) + H{-<j)i) - 2/7(0i) 

02 = Au;2 + if (03) + i^(-02) - if (-0l) - if (02) 

0, = Aio, + ff (0,+i) + H{-<p,) - ff (-0,-1) - H{<j),) 



!.^ = Acon + 2if (-0jv) - H{-<Pn-i) - HUm 
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(15) 



full model: g=.01 
Ti=6.0 



full model: g=.01 
Tl=6.0 
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Ti=6.0 



phase model: 
11=6.0 




I 



FIG. 2. Examples of traveling waves in both the Wang-Buszaki model (left) and the phase model 
(right) corresponding to two different values of the temperature dependent time constant r/. The 
phase model reproduces the dynamics of the full model. Furthermore, it is clear from the phase 
model, that as one increases the constant r/, the wavelength of the traveling waves decreases. 
Coupling strength in the full model is small: It must be on the order oi g = 0.01 ^S/cm? to 
ensure phase-locking. 

B. Traveling Waves In Chains Of Coupled Oscillators 



The wave solution for the travehng wave can be written: 



!)* + Awt, 



(16) 



where j is the index of the oscillator and (f)* is the phase shift between adjacent oscillators. 
Substituting this solution into the equations with non-reflecting boundary conditions, we 
see that dn] is a solution provided that H{(j)*) = H{—(f)*). Since our interaction function has 
both odd and even components, this statement is true only if 0* is the root of H{(j))odd- Thus, 
the stable phase-locked states for pair-wise symmetrically coupled oscillators determine the 
wavenumber, (p*. In a continuum limit, we can consider 6{x,t) to be a function of both 



11 



position and time. If we take the derivative of 6{x,t) with respect to t we obtain the 



expression 1 171]: 



de dedx ^ ,^^. 



which is equivalent to: 



U + (j)*Vg = 

.. = |i. (18) 

This is an expression for the phase velocity of the wave. Therefore, we see that we may 
identify the wavenumber of the system as: k = (p*. The stable phase- locked state between 
pairs of oscillators defines the wavenumber of a traveling wave. In the last section it was 
demonstrated that as we vary constants rj and Qk in the Wang-Buszaki model, the stable 
fixed point changes. This translates to a change in wavelength in a chain of neurons. Figure 
[2] shows a comparison between the phase model and full model for a variety of Qk and t]. 
The four panels on the left correspond to the full model. The panels on the right correspond 
to the phase model. The phase model quantitatively reproduces the dynamics of the full 
model. The phase model clearly demonstrates that the wavenumber increases for increasing 
7]. Coupling strength in the full model is small: it must be on the order of 0.01 to get 
agreement with the phase models. Stronger coupling results in only synchronous dynamics. 

C. Anti- Waves In Chains of Coupled Oscillators 

Anti-waves have been studied by Ermentrout and Kopell in two separate publications 
[l7||27|. The previous mechanisms for generating anti- waves rely on extremely long chains 
(essentially infinite) or chains with distal connections. 

Assuming an isotropic chain with no gradient in the natural frequencies, the intermediate 
stable phase- locked state defined by H{(f)*)odd = will generate traveling waves. If the fixed 
point 0* is identified as the wavenumber k, then the one kink anti-wave solution can be 
written: 

(j)j = kj + ut j < j* 

cP^ = -kj+cut j>j*. (19) 
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FIG. 3. Four examples of anti-waves in both rings and chains of oscillators computed with r/ = 6.0. 
A. Is a wave in a chain of Wang-Buszaki neurons with periodic boundary conditions. B. Is a 
wave in a chain of phase oscillators with periodic boundary conditions. C. Is a wave in a ring 
of Wang- Buszaki neurons with non-reflecting boundary conditions (full model). D. Is the phase 
model reduction of panel C: it is a wave in a chain of phase oscillators with non-reflecting boundary 
conditions. 

In these equations j* represents the position of the kink. By substituting the above expres- 
sion into[l5]we see again that the anti-wave is a solution provided that H{k) = H{—k) [k is 
the root of Hodd{k)) . Figure [3] demonstrates examples of anti-waves obtained in the Wang- 
Buszaki model compared with a phase model. Figures [3JA. and[3p are waves generated in the 
full model with non-periodic and non-reflecting boundary conditions, respectively. Figures 
EB and ED are the equivalent phase models. 

D. Obtaining Diflferent Wave solutions from Random Initial Conditions 

In order to see the variety of anti-waves, we will start chains of oscillators with random 
initial phases to estimate the basins of attraction [28i]. The interaction function we use f llip 
allows any number of shocks or kinks in the anti-wave, constrained only by the length of the 
chain. However, if we start the equations from random initial conditions, the probability of 
getting a certain number of shocks varies as we change the size of the Fourier components. 
The main point is that even if the odd portion of the interaction function allows for a multiple 
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FIG. 4. Probability of obtaining different solutions of equation (I15p for twenty oscillators as a 
function of the Fourier coefficients of the interaction function. A^ designates the number of shocks 
in an anti-wave. A^ = corresponds to a traveling wave. A. The probability distribution with 
H{(j)) = ai cos(0) + 6i sin((/)) + 0.75sin(2i?i>) as a function of ai with 6i = 1. B. The probability 
distribution with ai = 1 and varying bi. C. The probability distribution calculated with an 
interaction function: H{(j)) = —3 cos{(j)) — .92 cos{2(f)) + bi sin(0) — 0.75 sin(2i;^). D. The probability 
distribution calculated using the interaction function: H{(f)) = — 3cos(</>) + 6i sin(0) — 0.75sin(2(/)). 
For each parameter value the equations were solved from 10000 random initial conditions. 



shock solution, the probabihty of the system converging to this solution from random initial 
conditions may be extremely low and is determined by the magnitude of both the even and 
first odd Fourier modes. Figure H] shows the probability distributions of obtaining various 
anti-wave and traveling wave solutions as a function of the magnitudes of both the even 
and odd Fourier terms of the interaction function. Panel SJA. is a plot of the probability of 
obtaining either an A^-shock anti-wave solution or a traveling wave (0— shock) solution as 
a function of ai. From the plot, we see that the probability of obtaining a traveling wave 
solution approaches zero as ai — )■ and the probability distribution shifts towards A^ = 6. 
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FIG. 5. A. Compactons initiated from either side of the chain can shift the phase shock in the 
anti-wave back and forth. The interaction function used was: H{(j)) = cos{(j))+sva.{(j)) — 0.75 sin(2^) 
B. The shift in the kink for varying initial compacton amplitudes. 



As Oi increases towards 1, the probability of obtaining the traveling wave solution increases 
until it is the most probable state. Panel |1]B is a plot of the probability distribution as a 
function of hi with ai = 1. Panel HP shows a trend similar to|4]/\.. For 6i = the probability 
of obtaining a traveling wave solution from random initial conditions is close to 0. The 
solution with the maximal probability corresponds to an anti-wave with A^ = 9 shocks or 
kinks. As hi is increased towards 1 or decreased towards —1 the probability distribution 
shifts towards A^ = 0. That is, solutions with fewer kinks become more probable. Panels 
HP and HP demonstrate the effect of higher order even terms. Panel HP demonstrates that 
using the interaction function with two even terms {H{(j)) = cos(0) + hi sin(0) — |sin(20)) 
results in a probability distribution in which the traveling wave solution is the most probable 
solution for 6i — t- 1 and the 1 kink anti-wave solution is the most probable for hi — )■ —1. If 
we do not include this second order term, as demonstrated in panel HP the most probable 
solution corresponds to an (A^ = 4 shocks) anti-wave. The main point of this is that, even 
though only the odd Fourier modes determine the solutions to equations [13] and [151 the even 
Fourier modes can drastically affect the basin of attraction of those solutions. 
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E. Moving The Shock Position With Impulses 



Electric fish have been observed to produce complex anti-wave type patterns 27|] . What 
is more, the lead oscillator or (kink) in these anti-waves has been observed to be able 
to shift position. We present a mechanism which shifts the position of the kink/shock 
in the anti-wave without changing its shape. Over the last fifteen years, Pikovsky and 
Rosenau have written several papers studying waves in oscillator lattices with purely even 
coupling. They showed that such interaction functions can be derived from networks of 
Josephson junctions and Van-der Pol oscillators|36l][37[|38|. Waves in these networks take 
the form of solitary pulses which retain their shape. Pikovsky and Rosenau have coined 



these waves "compactons" [36| [37j |38| . Collisions between different compactons have been 
studied in detail. Nothing, however, has written about compactons colliding with a phase 
boundary (like the phase shock in the anti-wave). Compacton-like pulses are possible not 
only in systems with a purely even interaction function, but theycan be observed in systems 



with odd terms, provided the even component is large enough |38[|. The odd portion of the 
interaction function acts to dissipate the initial pulse, but a compacton-like wave may still 
travel large distances even with a substantial odd term. Compactons are possible in systems 
in which the interaction function generates anti-waves. For instance, we consider a system 
of 200 phase-difference equations with non-refiecting boundary conditions. The interaction 
function used is H{(j)) = cos(</)) + sin(</)) — .75 sin(20). The initial conditions are the one kink 
anti-wave solution with an extra pulse stimulus applied at the end. Thus we have the one 
kink/ shock solution: 

<t>M = k 3<- 

N 
0,(0) = ~k j> 



2 



(20) 



with the addition of a pulse: 



(pjiO) = k + -{l + cos{{j-xo)'K/a)) \j-xo\<a. (21) 

Here A is the amplitude, xq is the position of the pulse and a is the width of the pulse. In 
Figure [5] the pulse width used is a = 10 and \k\ = .84106. This pulse is the form used by 
Pikovsky et al. to generate compactons, but other initial conditions may work as well ^ 
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Figure [5p. shows multiple compactons traveling on top of an anti-wave and colliding with 
the shock located at A^ = 100. The anti-wave is composed of the stationary white and grey 
regions. The grey region of the plot corresponds to the solution = .841 whereas the white 
region corresponds to: (j) = 5.44. Upon collision, the shock shifts but retains its shape. In 
this manner, multiple pulses initiated at the ends of the chain may be used to shift the 
shock back and forth. Figure [5|3. is a plot of the shift in the kink as a function of the initial 
compacton amplitude. The maximum shift obtained is approximately 9 sites. li A > 1.5 
radians, larger amplitude pulses do not necessarily provide larger position shifts. If the pulse 
is too large, it will destroy the "perfect kink" solution. Thus, near the supercritical pitchfork 
bifurcation, not only can the shock of an anti-wave form anywhere along the the chain but 
precisely because of this property, it can be shifted around by an impulse (compacton). In 
this way, the additional even and odd Fourier terms produce a central pattern generator 
which is malleable. Perhaps this is a mechanism by which the hindbrain of a fish could send 
impulses to the rest of its spinal cord to modify the fish's swimming pattern. 

IV. STABILITY ANALYSIS 

In analyzing the stability of anti-wave solutions, one of the main questions we want to 
address is how the relative contributions of each Fourier mode contributes to the stability of 
the solution. Additionally, we want to examine the importance of other parameters in the 
model, such as the length of the chain and the position of the kink. In these chain models 
there are four basic wave solutions which are of interest. Two of the solutions correspond 
to a traveling wave in either direction (left to right or right to left) and two correspond 
to anti-waves. The first anti-wave solution describes a wave emanating from the center 
of the chain and propagating in both directions outward. The second anti-wave solution 
corresponds to two waves emanating from the edges and colliding in the center of the chain. 
Using our simplified interaction function we begin by analyzing the simplest (shortest) chain 
possible, the three neuron system. The three phase equations (122]) describing the neurons 
can be condensed to two by a change of variables. Linearizing these equations ([2SD about 
the anti-wave solution results in a 2x2 Jacobian. Thus, the problem is simple enough so that 
we can solve for the eigenvalues as a function of ai and subsequently show where and how 
the anti-wave solution loses stability. Once we have proved stability for this simple case, we 
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discuss longer chains and specifically, we analyze the effects of the magnitude of the even 
component on the stability of various anti-wave solutions. 

A. Stability Analysis Of The Three Oscillator System 

Our equations for the system with non-refiecting boundary conditions are: 

^'i = 2H{02 - 0,) 
62 = H{93 - 62) + H{9i - 62) 

4 = 2i7(^2-^3). (22) 

We then write them in terms of their phase differences: 

01 = 2/J(-02) - i^(-0i) - H{ct>2) 
02 = H{-<Pi) + i/(02) - 2iJ(0i). (23) 

Note that these equations are invariant under a reflection: 0i — )■ —02 

02 -^ -01 

Figure [6] is a plot of the nullclines of the system |23] for various values of ai using 61 = 1 
and 62 = —-75. There are two anti-wave solutions indicated by the boxes and two traveling 
wave solutions which are circled. One can see that when there is no even component the 
solutions to the system possesses perfect reflection symmetry. As soon as ai is nonzero, the 
system loses this symmetry. We want to analyze the stability of the anti-wave analytically. 
Choosing 0i = A; and 02 = —A;, we linearize our equations about this system and write down 
the Jacobian as follows: 

( H'{-k) -2H'{k)-H'{-k)\ ^ ^ 

Mo = \^ 24 

\^ -H'{-k) - 2H'{k) H'{-k) J 

Solving for the eigenvalues of this expression, we have 

Ai,2 = -2H'{k), -2H'{-k) - 2H\k). (25) 

As long as the derivative of these H functions evaluated at this solution is positive, then 
the solution will be stable. Substituting the simplified interaction function (|26|) and its 
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FIG. 6. Phase portrait of the two-dimensional system, ()23p where -ff(</>) = ai cos((/>) + sine/; — 
0.75 sin 20. The 0i (^2) nuhchne is in red(green). A. The nullcHnes and the flow for a\ = 0. The 
fixed points corresponding to anti-waves are enclosed with boxes. The fixed points corresponding 
to traveling waves are circled. B. The nullclines and the fiow for ai = 1.2. C. Bifurcation diagram 
computed with AUTO: The anti-wave solution loses stability near oi = 1.118. In this plot, the 
stable solution is represented by a solid line whereas the unstable solution is represented by the 
dotted lines. 

derivative into equation [23 results in the eigenvalue expressions: equation |271 

H{(f)) = hi sin(0) + 62 sin(20) + ai cos(0) 
H\(l)) = hi cos(0) + 262 cos(20) - ai sin(0) 

(26) 



Ai = — 2(61 cos(/c) + 262 cos(2/c) — ai sin(x)) 
A2 = -4(61 cos(A;) + 262 cos(2A;)) 



(27) 



At the critical value of the parameter acriucai = l-l? the first eigenvalue vanishes. For 
larger values of ai, the fixed point becomes unstable. The mirror symmetry of the equations, 
along with the bifurcation dia gram in Figure [6l suggest that this is a subcritcal pitchfork 
bifurcation. This is verified in |39| where we explicity calculate the normal form equations. 
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B. Stability Analysis For The Traveling Wave 

Proving the stability of the travehng wave (with no kinks) is relatively straightforward. 
We begin by writing down the equations for the "jth" oscillator. Since the system is nearest 
neighbor coupling, the general expression is: 

<jy, = H{<P^^,) + H{-<P,) - i/(-</.,-i) - H{<j),). (28) 

Defining a = if (0*), b = H' {—</)*) and linearizing the equations about the wave solution, 
the equations for the phase are simply a discretized version of Laplace's equation: 

a(/)j+i — (a + b)(f)j + b(t>j-i = X(f)j. (29) 

We may solve for the equations by assuming a general solution: 

(pj = Axi + Bxi, (30) 

and invoking the boundary conditions: (po = (f>i and (pN-i = (pN .Plugging our solution [30] 
into Eni we may solve for the eigenvalues of the system: 

3fJA = 2v^cos-^^-(a + 6). (31) 

Therefore the wave solution will always be stable provided that 

2Vba <a + b (32) 

Or alternately that a,b > 0. Stability is lost for a,b < 0. 

C. Stability for The Anti-Wave under More General Conditions 

Stability analysis of the anti-wave solutions may be proven using a combination of the 
Gershgorin Circle Theorem and numerically computing the eigenvalues of the linearized 
equations. In these examples, we use non-reflecting boundary conditions. The argument will 
be identical for periodic boundary conditions. One starts by assuming a 1 shock solution 
and linearizing the phase the equations about this solution. The Jacobian will have matrix 
elements of the form: 

a,, = -{H'{<p,) + H'{ct>,)) a,,,+i = i/'(0,+i) a,,_i = if'(0,_i) (33) 

20 



All other matrix elements are zero. Once more, define: a = H'{(j)*), b = H'{—(j)*). If / 
denotes the location of the kink, then the matrix elements corresponding to the kink are: 

ai,i = -{a + b) ai^i+i = b ai^i^i = b (34) 

or, 

0-1,1 = -{a + b) ai^i+i = a ai^i-i = a (35) 



depending on the orientation of the shock. The Gershgorin Circle Theorem J40| tells us 
that all of the eigenvalues of the matrix lay in the union of disks centered at the diagonal 
elements of the matrix with radii less than the absolute value of the sum of the row entries 
not including the diagonal terms. If the Jacobian of the anti-wave is described by equation 
|3l]all eigenvalues will lie in the union of three disks: one centered at — (a + b) with radius 
a + b, one centered at —2a + b (corresponding to the ends of the chain) with radius a and 
one centered at — (a + b) with radius 26. If we assume that the even term ai is positive and 
increasing, then H'[—(f)) > H'{(j)) 

\b\ > \a\, 
|26|>|a + 6|. (36) 

Thus the disks will extend beyond the origin, and we will not be able to say anything about 
stability using this theorem. On the other hand, if the solution is a shock oriented in the 
opposite direction, such as in equation |35l the disk corresponding to the equation at the 
shock is centered at —{a + b) and extends out with radius \2a\. The disk will always lie in the 
left half of the imaginary plane for ai positive. Therefore, this solution will always be stable. 
In the former case, in which the fractured wave is not necessarily stable, one must apply 
numerical methods to explicitly calculate the eigenvalues of the Jacobian. We are interested 
how variables such as the position of the shock and the length of the chain effect the stability 
of the solution as we vary the magnitude of the first even Fourier mode. Again, we use the 
interaction function: H{(f)) = bi sm{(f)) + 62 sin(20) + ai cos{(f)), where 61 = 1 and 62 = —.75. 
We start by examining a fifty-one oscillator chain (described by fifty equations) in which we 
move the position of the shock. Figure [TK shows the critical value of ai at which the shock 
solution loses stability as a function of shock position where the position varies between site 
4 and site 46. The parameter a critical is determined by calculating the eigenvalues of the 
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FIG. 7. A. Critical value of a\ as a function of shock position. B. Eigenvalue with the maximal 
real part as a function of a\ for different shock positions. N denotes the location of the shock. C. 
An interaction function with a more negative 62 can also possesses a larger a\ before the solution 
becomes unstable. 

Jacobian for various values of a^ and determining when the eigenvalue with the maximal 
real part becomes positive. Panel [7B plots the real part of the eigenvalue with maximal 
real part as a function of ax. as the position shifts, the eigenvalues lose stability at different 
values of a\. The type of bifurcation by which they lose stability changes as well. For 
shocks located at even numbered sites, the system loses stability in a Hopf bifurcation as a 
complex conjugate pair of eigenvalues crosses the origin simultaneously. For shocks located 
at odd sites the system apparently loses stability in a sub-critical pitchfork bifurcation (This 
has not been rigorously proven). Figure [TP is a plot of acriucai for the eigenvalues of the 
fifty-one oscillator Jacobian linearized about solutions corresponding to varying values of 62- 
The plot clearly demonstrates that shock solutions corresponding to a larger value of 62 can 
support a larger even component before becoming unstable. Figure [8] is a plot of acriucai vs. 
the number of phase difference equations (number of oscillators in a chain). The shock is 
always located centrally: e.g. for 101 oscillators, the shock is located at A^ = 50. For an 
odd number of phase difference equations, the shock is position is obtained by dividing the 
number of equations in half and rounding up. Figure [8] shows phenomena similar to Figure 
lOA. That is, a solution will lose stability for different a\ dependent on the position of the 
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phase difference equation number 

FIG. 8. Critical value of ai as a function of chain length where N represents the number of phase 
difference equations (J)n- Depending on the length of the chain and position of the shock, the 
solution may lose stability in either a Hopf bifurcation or what is believed to be a subcritical 
pitchfork bifurcation. The shock is always located centrally: (e.g. for 101 oscillators, the shock 
is located at N=50. For an odd number of phase difference equations, the shock is position is 
obtained by dividing the number of equations by half and rounding up.) 

shock in the solution as well as the number of oscillators in the chain. In the case that the 
shock is perfectly centered, we expect the solution to lose stability in a sub-critical pitchfork. 
Regardless of the position of the shock, even for relatively short chains the anti-wave solution 
will be stable for a relatively large even component, whose magnitude is at least as large as 
the first odd Fourier mode. The manner in which the solution loses stability and the size 
of the even component it can support depend on both the position of the kink/shock and 
the length of the chain. The Gershgorin circle theorem tells us that one of the anti-wave 
solutions will always be stable, no matter how long the chain. 

V. PATTERN FORMATION IN TWO-DIMENSIONAL ARRAYS 



Anti-wave patterns can be observed in two-dimensional networks as well. Using nearest 
neighbor coupling, the differential equation for a single oscillator is: 

,, = H{9x+l,y — 9x,y) + H{9x,y+l — Ox,y) + H{9x,y-l — 9x,y) + H{9x-l,y — 9x,y), (37) 

where x, y are discrete indices describing the location of an oscillator. These indices run 
from 1 to A^ where N"^ is the number of differential equations in the array. As in the one 
dimensional case, we use non-reflecting boundary conditions. Examples of these patterns 
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FIG. 9. Examples of two-dimensional patterns. The horizontal and vertical axis are the oscil- 
lator indices. A. Snapshot of a slowly varying wave pattern obtained from compacton-like initial 
conditions: a 2-dimensional pulse is initiated in the upper left-hand corner of the array on a 
synchronous background. The interaction function used is H{(f)) = sin((/)) — cos(i;^) — .75sin(2i?i)) 
B. Stable stationary fractured pattern obtained from random initial conditions and interaction 
function:— 2 cos(i;^) — 0.518 sin((/)) — 1.31sin(2(/>) — .933sin(30) C. Stationary anti-wave generated 
with H{(f)) = cos{(j)) + sin((/)) — .75sin(2ij!)). D. Stationary traveling wave generated with the same 
interaction function as C. 

and the interaction functions used to generate them are are illustrated in Figure [H Patterns 
|9j3j9p and[9p are stationary and numerically stable. Pattern [9l3 is generated from random 
initial conditions. Patterns HP and |9p are stable but have a very small basin of attraction. 
The initial conditions must be almost identical to the actual anti-wave or traveling wave 
solution. 

Analogous to the one-dimensional case, a plane wave may be represented by the solution: 



Oxy = Kx + kyy + cot. 
Whereas the shock solution may be written: 



(38) 



^xy '^X'^ ~i~ '^yU 



9. 



xy 



-KxX ~r kylj 



X < X 
X > X . 



(39) 



In these equations x* is the location of the shock. A shock may be thought of as a discon- 
tinuous boundary between the two traveling waves. The fractured pattern in Figure |9}3 is 
composed of many small waves of varying kx and ky which form shocks. Fractured patterns 
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are generated with random initial conditions where the phases are chosen between and 



27r|4l| |42| 



VI. CONCLUSION 

Intermediate stable phase-locked states occur in a variety of neuron models and are 
capable of producing a pattern of wave-like activity known as an anti-wave. This type of 



activity was first observed in the dogfish spinal cord by Grillner 32| and it may occur in other 



networks as well[33|][l2|. The mechanism we use for generating these waves may have an 
experimental basis: Peinado et al. have shown that modulating the potassium conductance 
of gap junction coupled neurons in the neonatal mouse neocortex leads to wave behavior 



13| . The mechanism for generating traveling waves and anti- waves is more flexible than past 
models, because it allows one to modulate the wavelength of the wave by adjusting properties 
such as the potassium conductance or temperature dependent time constant. It allows for 
a huge variety of phase-locked patterns in a chain or array of neurons. In the case of anti- 
wave generation, it does not require distal connections and the shock can form anywhere 
along the chain. The shock position can also be moved by colliding phase compactons. This 
may provide a mechanism by which an animal such as the dogfish could switch between 
different swimming patterns. More generally, near the bifurcation from in-phase synchrony, 
the interaction functions of the phase model have large higher order Fourier modes. This is in 
contrast to many oscillator models that were used in the past, which only include the lowest 
odd Fourier mode. The odd modes determine the value of the stable-phase locked solution 
but the even modes affect the stability of the anti-wave solution. Varying the relative even 
component affects the basin of attraction of a particular solution, or the probability that 
the phases will converge to a particular anti-wave solution from random initial conditions. 
Finally, we note that this behavior is not only relevant to chains but to two dimensional 
arrays as well. 

Appendix A: The Wang-Buszaki Model 

The parameters used in the Wang-Buszaki model are: Vgyn = — 60.5mV, ql = 0.1/iS, 
vl = — 65mV, gNa = 35 fiS/cm^, Vatq = 55mV, gx = ^ fiS/cm?, Vk = — 90mV, Ojo = 4 
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Ti = ISms and io = .63 nA/cm? . The nonlinearites mentioned in Section fll Al are given by 
the following expressions: 



a^{V) = 0.l(y + 35.0)/(1.0 - exp(-(y + 35.0)/10.0)) 

/3^(V) = 4.0exp(-(V + 60.0)/18.0) 

M^V) = aUV)/{am{V) + P^V)) 

ah{V) = 0.07exp(-(\/ + 58.0)/20.0) 

Ph{V) = 1.0/(1.0 + exp(-(\/ + 28.0)/10.0)) 

H^{V) = ah{V)/{ah{V) + f3f,{V)) 

a^(V) = 0.01(\/ + 34.0)/(1.0 - exp(-(V" + 34.0)/10.00)) 

(3n{V) = 0.125 exp(-(\/ + 44.0)/80.0) 

N^V) = an{V)/{an{V) + I3^{V)) 
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[41] All equations in this section were integrated with Euler's method. The step size used is .01. 

[42] The Fourier terms used may or may not be calculated from the full model. For instance, an 
interaction function with an artificially inflated second odd mode will produce an interesting 
fractured pattern, whereas many interaction functions computed from the full model will not. 



[43] D. Somers and N. Kopeh, Physica D: Nonlinear Phenomena 89, 169 (1995), ISSN 0167-2789, 



|http : //www. sciencedirect . com/science/article/pii/01 67278995001980 

[44] I. Gradshteyn and I. Ryzhik, Tables of Integrals, Series and Products (Academic Press, 2000) 
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